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Confined suspensions of active particles show peculiar dynamics characterized by wall ac¬ 
cumulation, as well as upstream swimming, centerline depletion and shear-trapping when 
a pressure-driven flow is imposed. We use theory and numerical simulations to investigate 
the effects of confinement and non-uniform shear on the dynamics of a dilute suspen¬ 
sion of Brownian active swimmers by incorporating a detailed treatment of boundary 
conditions within a simple kinetic model where the configuration of the suspension is de¬ 
scribed using a conservation equation for the probability distribution function of particle 
positions and orientations, and where particle-particle and particle-wall hydrodynamic 
interactions are neglected. Based on this model, we first investigate the effects of con¬ 
finement in the absence of flow, in which case the dynamics is governed by a swimming 
Peclet number, or ratio of the persistence length of particle trajectories over the channel 
width, and a second swimmer-specific parameter whose inverse measures the strength 
of propulsion. In the limit of weak and strong propulsion, asymptotic expressions for 
the full distribution function are derived. For finite propulsion, analytical expressions for 
the concentration and polarization profiles are also obtained using a truncated moment 
expansion of the distribution function. In agreement with experimental observations, the 
existence of a concentration/polarization boundary layer in wide channels is reported 
and characterized, suggesting that wall accumulation in active suspensions is primarily a 
kinematic effect which does not require hydrodynamic interactions. Next, we show that 
application of a pressure-driven Poiseuille flow leads to net upstream swimming of the 
particles relative to the flow, and an analytical expression for the mean upstream velocity 
is derived in the weak flow limit. In stronger imposed flows, we also predict the forma¬ 
tion of a depletion layer near the channel centerline, due to cross-streamline migration 
of the swimming particles towards high-shear regions where they become trapped, and 
an asymptotic analysis in the strong flow limit is used to obtain a scale for the deple¬ 
tion layer thickness and to rationalize the non-monotonic dependence of the intensity 
of depletion upon flow rate. Our theoretical predictions are all shown to be in excellent 
agreement with finite-volume numerical simulations of the kinetic model, and are also 
supported by recent experiments on bacterial suspensions in microfluidic devices. 
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1. Introduction 

The interaction of active self-propelled particles with rigid boundaries under confine¬ 
ment plays a central role in many biological processes. Spermatozoa are well known to 
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accumulate at rigid boundaries (Rothschild 1963; Woolley 2003), with complex implica¬ 
tions for their transport in the female tract during mammalian reproduction (Suarez & 
Pacey 2006; Denissenko et al. 2012; Kantsler et al. 2014). The aggregation of bacteria 
near surfaces and their interaction with external flows in confinement has a strong effect 
on their ability to adhere and form biofilms (Rusconi et al. 2010; Lecuyer et al. 2011; Kim 
et al. 2014). It also impacts their interactions with the gastrointestinal wall during diges¬ 
tion, with consequences for various pathologies (Lu & Walker 2001; Cellia et al. 2009). 
Confinement has also been shown to affect cell-cell interactions and collective motion in 
dense sperm and bacterial suspensions and can also result in spontaneous unidirectional 
flows (Riedel et al. 2005; Wioland et al. 2013; Lushi et al. 2014). In engineering, the 
ability to concentrate or separate bacteria by controlling their motions in microfluidic 
devices with complex geometries has been demonstrated (Galajda et al. 2007; Hulme 
et al. 2008; Lambert et al. 2010; Kaiser et al. 2012; Altshuler et al. 2013), as well as 
the ability to harness bacterial swimming power to actuate gears (Sokolov et al. 2010; 
Di Leonardo et al. 2010) or transport cargo (Koumakis et al. 2013; Kaiser et al. 2014). 
Particle-wall interactions are also critical in systems involving synthetic microswimmers 
(Gibbs et al. 2011; Takagi et al. 2013, 2014), as these inherently reside near surfaces due 
to sedimentation. 

The prominent feature of confined active suspensions is the tendency of swimming 
particles to accumulate near boundaries. This was first brought to light by Rothschild 
(1963), who measured the concentration of swimming bull spermatozoa in a glass cham¬ 
ber and reported a nonuniform distribution across the channel with a strong spike in 
concentration near the walls. Berke et al. (2008) repeated the same experiment using 
suspensions of Escherichia coli in microchannels and also observed an accumulation of 
bacteria at the channel walls. They further reported the tendency of bacteria to align 
parallel to the boundaries, which led them to consider wall hydrodynamic interactions 
due to the force dipole exerted on the fluid by the self-propelled particles as a poten¬ 
tial mechanism for migration. Hydrodynamic interactions are indeed known to have an 
impact on the trajectories of swimming particles near no-slip walls (Lauga et al. 2006; 
Spagnolie & Lauga 2012), and have been shown to lead to attraction of sperm cells 
towards walls (Fauci & McDonald 1995). Li & Tang (2009) and Li et al. (2011) also 
observed wall accumulation in suspensions of Caulobactor crescentus but presented an 
alternate mechanism based purely on kinematics that explains accumulation as a result 
of the collisions of the bacteria with the wall, leading to their reorientation parallel to 
the surface. The possibility of a non-hydrodynamic mechanism for wall accumulation is 
indeed supported by various simulations that neglected wall hydrodynamic interactions 
(Gostanzo et al. 2012; Elgeti & Gompper 2013), suggesting that such interactions in fact 
only play a secondary role in this process. 

Several other interesting effects have also been reported when an external flow is ap¬ 
plied on the suspension. One such effect is the propensity of motile particles to swim 
upstream in a pressure-driven flow. This was noted for instance by Hill et al. (2007), 
who tracked the trajectories of Escherichia coli in a shear flow near a rigid surface in a 
microfluidic channel, and proposed a complex mechanism for upstream swimming based 
on the chirality of the flagellar bundles and on hydrodynamic interactions. Such inter¬ 
actions were characterized more precisely by Kaya & Koser (2009), who demonstrated 
that the Escherichia coli cells undergo modified Jeffery’s orbits (Jeffery 1922) near the 
walls and suggested that this detail is crucial in understanding the upstream migration. 
A clearer picture of this phenomenon emerged in yet more recent work by Kaya & Koser 
(2012), who systematically analyzed Escherichia coli motility near a surface as a func- 
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tion of the local shear rate. At low shear rates, circular trajectories were observed due 
to the chirality of the cells, as previously explained by Lauga et al. (2006). At higher 
shear rates, positive rheotaxis was reported and accompanied by rapid and continuous 
upstream motility. This directional swimming was explained as a result of the combined 
effects of surface hydrodynamic interactions, which were thought to cause the swimming 
cells to dip towards the walls, and of reorientation by the shear flow, which aligns the cells 
against the flow. Upstream motility was also recently discussed by Kantsler et al. (2014) 
in the case of mammalian spermatozoa, where the combination of shear alignment, wall 
steric interactions and cell chirality was shown to lead to steady spiraling trajectories in 
cylindrical capillaries. 

While most experimental studies under confinement have focused on near-wall aggre¬ 
gation and swimming dynamics, the behavior of self-propelled micro-organisms under 
flow in the bulk of the channels is also of interest. In recent work, Rusconi et al. (2014) 
analyzed the effects of a Poiseuille flow on the trajectories and distributions of motile 
Bacilus subtilis cells, with focus on the central portion of the channel. In sufficiently 
strong flow, they reported the formation of a depletion layer in the central low-shear re¬ 
gion of the channel, accompanied by cell trapping in the high-shear regions surrounding 
the depletion. This trapping was attributed to the strong alignment of the swimming cells 
with the flow under high shear, which hinders their ability to swim across streamlines. 
Quite curiously, they reported that maximum depletion is achieved at a critical imposed 
shear rate of approximately 10 s“^, above which both trapping and depletion become 
weaker. A simple Langevin model capturing the effects of self-propulsion, shear rotation, 
and diffusion was also proposed to explain these observations, and was able to reproduce 
the salient features of the experiments. 

Models and simulations explaining the mechanisms leading to these rich dynamics 
have been relatively scarce. Direct numerical simulations of hydrodynamically interact¬ 
ing swimming particles confined to a gap between two plates were first performed by 
Hernandez-Ortiz et al. (2005, 2009) using a simple dumbbell model, and indeed cap¬ 
tured a strong particle accumulation at the boundaries in dilute systems. As the mean 
swimmer density was increased, collective motion and mixing due to particle-particle 
hydrodynamic interactions led to a decrease in the concentration near the walls. Accu¬ 
mulation was also observed in simulations of self-propelled spheres by Elgeti & Gomp- 
per (2013), who entirely neglected hydrodynamic interactions. This study, as mentioned 
above, suggests that wall hydrodynamic interactions are not required to explain migra¬ 
tion, and neither is shape anisotropy. Rather, the simple combination of cell swimming, 
steric exclusion by the walls, and diffusive processes is sufficient to capture accumulation, 
and Elgeti & Gompper (2013) also proposed a simple Eokker-Planck description of the 
suspension that shares similarities with the present work and was able to explain their 
results. A similar continuum model was also proposed by Lee (2013), who derived ana¬ 
lytical expressions for the ratio of particles in the bulk vs near-wall region in the limits 
of weak and strong rotational diffusion. Very recently, Li & Ardekani (2014) performed 
direct numerical simulations of confined suspensions of spherical squirmers that propel 
via an imposed slip velocity, and reported strong accumulation at the boundaries irre¬ 
spective of the details of propulsion. They also noted the tendency of particles to align 
normal to the wall in the near-wall region. 

The effects of an external flow have also been addressed using discrete particle mod¬ 
els and simulations. The dynamics of isolated deterministic microswimmers in Poiseuille 
flow were studied in detail by Zdttl & Stark (2012, 2013), who found that such swimmers 
perform either an upstream-oriented periodic swinging motion or a periodic tumbling 
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motion depending on their location in the channel. Suspensions of interacting swimmers 
in pressure-driven flow have also been simulated, notably by Nash et al. (2010) and 
Costanzo et al. (2012), who both observed aggregation at the walls together with up¬ 
stream swimming as a result of the rotation of the particles by the flow. More recently, 
Chilukuri et al. (2014) extended the simulation method of Hernandez-Ortiz et al. (2009) 
to account for a Poiseuille flow. Similar trends as reported earlier were observed, includ¬ 
ing wall accumulation and upstream swimming, as well as the reduction of accumulation 
with increasing flow rate. In addition, they also reported the formation of a depletion 
layer near the channel centerline in strong flows, in agreement with the microfluidic ex¬ 
periments of Rusconi et al. (2014). Simple scalings for the dependence of this depletion 
with shear rate, swimming speed and channel width were also proposed. 

While these various numerical simulations have been able to reproduce the relevant 
features of previous experiments, a clear unified theoretical model capable of capturing 
and explaining all of the above effects based on conservation laws and microscopic swim¬ 
mer dynamics is still lacking. In unconfined systems, much progress has been made over 
the last decade in the description of the behavior of active suspensions using continuum 
kinetic theories (Saintillan & Shelley 2013; Marchetti et al. 2013; Subramanian & Koch 
2009). One such class of models, introduced by Saintillan & Shelley (2008o,&) to explain 
the emergence of collective motion in semi-dilute suspensions, is based on a conservation 
equation for the distribution function 'P{x,p,t) of particle positions and orientations, in 
which fluxes arise due to self-propulsion, advection and rotation by the background fluid 
flow, as well as diffusive processes. When coupled to a model for the fluid flow (whether 
externally imposed or driven by the swimmers themselves), this conservation equation 
can be linearized for the purpose of a stability analysis or integrated in time to inves¬ 
tigate nonlinear dynamics. This approach, which also relates to other models developed 
in the context of active liquid crystals (Baskaran & Marchetti 2009; Marchetti et al. 
2013; Forest et al. 2013), has been very successful at elucidating the mechanisms leading 
to collective motion at a suspension level. However, attempts to apply such continuum 
kinetic theories to confined suspensions have been few and far between, in part due to 
the complexity of the boundary conditions that need to be enforced on the distribution 
function. 

In this paper, we present a simple continuum theory for the dynamics and transport 
of a dilute suspension of Brownian active swimmers in a pressure-driven channel flow 
between two parallel flat plates. To focus on the effects of steric confinement and its 
interaction with the flow, we neglect particle-particle and particle-wall hydrodynamic 
interactions entirely but incorporate a detailed treatment of the boundary conditions 
for the distribution function. As we show below, our theory is able to capture all the 
different regimes discussed above, including wall accumulation in the absence of flow, 
and upstream swimming, depletion at the centerline and trapping in high-shear regions 
when a flow is applied. We introduce the governing equations, boundary conditions and 
nondimensionalization in §2, where we also derive a simpler approximate model based on 
moment equations. The equilibrium distributions in the absence of flow are obtained in 
§3, where wall accumulation is seen to be accompanied by a net polarization of the par¬ 
ticle distribution near the boundaries, and where a very simple expression is derived for 
the concentration profile across the channel in terms of the parameters of the problem. 
The effects of an external Poiseuille flow are discussed in §4, where a numerical solu¬ 
tion of the governing equations captures upstream swimming and shear trapping in the 
relevant parameter ranges, and where both effects are also explained theoretically using 
asymptotic analyses in the weak and strong flow regimes. We summarize our results in 
§5 and discuss them in the light of the recent literature in the field. 
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Figure 1. Problem definition: a dilute suspension of slender active particles with positions 
X = {x,y,z) and orientations p = (sin 0 cos c;/), sin 6^ sin 0, cos 0) is confined between two parallel 
flat plates (2 = ±_ff) and subject to an imposed pressure-driven parabolic flow. 

2. Governing eqnations 

2.1. Problem definition and kinetic model 

We analyze the dynamics in a dilute suspension of self-propelled slender particles confined 
between two parallel flat plates and placed in an externally imposed pressure-driven flow 
as illustrated in figure 1. The channel half-width is denoted by iJ, and is assumed to 
be much greater than the characteristic length L of the particles {H/L ^ 1), so that 
the finite size of the particles can be neglected. The external flow follows the parabolic 
Poiseuille profile 

U{x) = U{z)y = Um[l-{z/Hf]y, (2.1) 

with maximum velocity Um at the centerline (2 = 0). The shear rate varies linearly with 
position 2 across the channel: 

= ( 2 . 2 ) 

where fiw = “^Um/H is the maximum absolute shear rate attained at the walls (2 = ±H). 

Following previous models for active suspensions (Saintillan & Shelley 2008a,6), the 
configuration of the active particles is captured by the probability distribution function 
'T{x,p,t) of finding a particle at position x = {x,y,z) with orientation p = (sin 0 cos 
sin0sin(/), COS0) at time t, where p also defines the direction of swimming. Conservation 
of particles is expressed by the Smoluchowski equation (Doi & Edwards 1986) 

d<P 

— + V,-{xT) + Vp-{pT)=0, (2.3) 

where the translational flux velocity x captures self-propulsion with constant velocity 14 
in the direction of p, advection by the imposed flow, and center-of-mass diffusion with 
isotropic and constant diffusivity dp. 

x = V,p + U{z)-dtV^\nT. (2.4) 

Particle rotations are captured by the angular flux velocity p, which includes contribu¬ 
tions from the imposed flow via Jeffery’s equation (Jeffery 1922; Bretherton 1962), and 















6 


B. Ezhilan and D. Saintillan 


from rotational diffusion with diffusivity dr'. 

p = S{z){z‘p){l -pp)-y- drVplnW. (2.5) 

We have assumed that the particles have a high aspect ratio, a good approximation for 
common motile bacteria as well as many self-propelled catalytic micro-rods. Particle- 
particle hydrodynamic interactions have also been neglected based on the assumption of 
infinite dilution; such interactions could otherwise be included via an additional distur¬ 
bance velocity in the expressions for x and p (Saintillan & Shelley 20086). As a result, 
we expect the distribution of particles to be uniform along the x and y directions, and at 
steady state the Smoluchowski equation (2.3) for E{x,p,t) = 'P{z,p) then simplifies to 

VsCosO— - + S{z)Vp • [cosd(/ -pp) • y<6'] = dr^l^. (2.6) 

This equation simply expresses the balance of self-propulsion, translational diffusion, 
particle alignment by the imposed flow, and rotational diffusion. 

In this work, we treat the translational and rotational diffusivities dt and dr as inde¬ 
pendent constants, which could result from either Brownian motion or various athermal 
sources of noise (Drescher et al. 2011; Garcia et al. 2011). The athermal contribution 
to diffusion may arise due to tumbling or other fluctuations in the swimming actuation 
of motile micro-organisms, or from fluctuations in the chemical actuation mechanism of 
catalytic particles. In many active suspensions, such athermal fluctuations are in fact the 
dominant source of diffusion. 


2.2. Boundary conditions 

In the continuum limit, the impenetrability of the channel walls is captured by prescribing 
that the normal component of the translational flux be zero at both walls: 

z • X = 0 at z = ±H. (2-7) 

Inserting equation (2.4) for the translational flux, this leads to a Robin boundary condi¬ 
tion for the probability distribution function: 

dE 

dt^^ = VaCos6'B at z = ±H, (2.8) 

oz 

expressing the balance of translational diffusion and self-propulsion in the wall-normal 
direction. Equation (2.8) implies that particles pointing towards a wall (cos 6* > 0 for 
the top wall at z = -\-H) incur a positive wall-normal gradient {dE/dz > 0), whereas 
particles pointing away from the wall {cos 6 < 0) incur a negative gradient. This suggests 
that sorting of orientations should occur and lead to a net polarization towards the walls, 
accompanied by near-wall accumulation. These effects will indeed be confirmed in §3. It 
is important to note that the boundary condition (2.8) requires that the wall-normal 
swimming flux be balanced by a diffusive flux. In the complete absence of translational 
diffusion {dt = 0), the swimming flux can no longer be balanced at the wall: this singular 
limit, which is ill-posed in our mean-field theory, will not be addressed here. Note also 
that the balance of the wall-normal fluxes hints at a length scale of = dt/Vg for wall 
accumulation, as we demonstrate more quantitatively below. 

Other types of boundary conditions have been considered in previous works. In par¬ 
ticular, several studies have implemented the condition 


z • xEdp = Q at z = ±i7. 


(2.9) 



7 


Transport of a dilute active suspension 

where 12 denotes the unit sphere of orientation. Equation (2.9) captures the zeroth ori¬ 
entational moment of (2.8) and is easily implemented numerically using a reflection 
condition on the distribution function. It was first used by Bearon et al. (2011) in a 
two-dimensional model of suspensions of gyrotactic swimmers constrained to a planar 
domain. Ezhilan et al. (2012) also imposed equation (2.9) in the case of a chemotactic 
active suspension confined to a thin liquid film, where the primary mechanism for ac¬ 
cumulation was chemotaxis as opposed to kinematics. In the absence of external fields, 
however, this boundary condition allows for a uniform isotropic solution throughout the 
channel and is therefore unable to capture near-wall accumulation or upstream swimming 
when a flow is imposed (see Appendix A for more details). Kasyap & Koch (2014) also 
considered chemotactic active suspensions in thin films but used a position/orientation 
decoupling approximation for the probability distribution function T {x,p,t), allowing 
them to derive a boundary condition for the number density field expressing the balance 
of the chemotactic and diffusive fluxes at the boundaries. To our knowledge, the only 
previously reported use of the boundary condition (2.8) for a confined active suspension 
was in the work of Elgeti & Gompper (2013), whose analysis was restricted to equilib¬ 
rium distributions in the absence of flow and in the limits of narrow channels or weak 
propulsion. 

Finally, it should be kept in mind that the simple boundary condition (2.8) neglects the 
finite size of the particles and is therefore inaccurate very close to the walls, where steric 
exclusion prohibits certain particle configurations and should lead to a depletion layer 
as observed in experiments (Takagi et al. 2014). The implications of steric exclusion are 
discussed further in Appendix B, where a more detailed boundary condition is derived and 
enforced on the hypersurface separating allowed from forbidden configurations (Nitsche & 
Brenner 1990; Schiek & Shaqfeh 1995; Krochak et al. 2010). As we show there, the effects 
of steric exclusion are weak in wide channels {H/L 3> 1) such as the ones considered in 
this work. 


2.3. Dimensional analysis and scaling 

Dimensional analysis of the governing equations reveals three dimensionless groups: 


PCs = 


Vs 


2drH' 


Y-t 'Iw A dt dr 

= ‘'‘vF- 


( 2 . 10 ) 


The first parameter Pe^, or swimming Peclet number, can be interpreted as the ratio of 
the characteristic timescale for a particle to lose memory of its orientation due to rota¬ 
tional diffusion over the time it takes it to swim across the channel width. Equivalently, 
it is also the ratio of the persistence length of particle trajectories {ip = Vg/dr) over 
the channel width (2iJ). The second parameter Pe/, or flow Peclet number, compares 
the same diffusive timescale to the characteristic time for a particle to align under the 
imposed velocity gradient. The third parameter A relates the translational and rotational 
diffusivities to the swimming speed and is a fixed constant for a given particle type. It 
can be interpreted as an inverse measure of the strength of propulsion of a swimmer with 
respect to fluctuations, and the limits of A —)■ 0 and A —> oo describe the strong and weak 
propulsion cases, respectively. When A is held constant, Pe^ also reduces to an inverse 
measure of confinement, with Pcg —>■ 0 and PCg —)■ oo describing the limits of weak and 
strong confinement, respectively. 

In the following, we nondimensionalize the governing equations using the characteristic 
time, length and velocity scales 


tc dj. , 


£c = H, Vc = Hdr, 


( 2 . 11 ) 
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and also normalize the distribution function W by the mean number density n defined as 

n= ^ J ^ J^<I'{z,p)dpdz. (2.12) 

After nondimensionalization, the conservation equation (2.6) becomes 

Pescosd^ - 2APe1^^ + ^ S{z)'Vp • [cos6»(/ -pp) -yE] = (2.13) 

where the dimensionless shear rate profile is simply S{z) = —z. The boundary condition 
(2.8) also becomes 

1 

^—= —T- COS0E at z = ±1. (2-14) 

dz 2APes 

Note that the choice of H for the characteristic length scale is convenient as it sets the 
positions of the boundaries to z = ±1 in the dimensionless system. However, we will see 
below that alternate length scales are more judiciously chosen in certain limits due to 
the presence of boundary layers. 

2.4. Orientational moment equations 

Equation (2.13), together with boundary condition (2.14), cannot be solved analytically 
in general. While a numerical solution is possible as we show below, analytical progress 
can still be made in terms of orientational moments of the distribution function (Saintillan 
& Shelley 2013). More precisely, we introduce the zeroth, first, and second moments of 
'P{z,p) as 

c(z) = (l), m{z) = {p), D{z) = {pp-1/3), (2.15) 

where the brackets (•) denote the orientational average 

(hip)) = [ h(p)tf'(z,p)dp. (2.16) 

Jn 

The zeroth moment c(z) corresponds to the local concentration of particles. The next two 
moments are directly related to the polarization vector P{z) and to the nematic order 
parameter tensor Q{z) commonly used in the description of liquid-crystalline systems 
(Marchetti et al. 2013) as 

Tn(z) = c(z)P(z), D{z) = c(z)0(z). (2-17) 

Knowledge of these as well as higher moments also allows one to recover the full distri¬ 
bution function as 

13 15 

p) = ^ c(z) + ^P‘ + —pp: D{z) + ..., (2.18) 

which can also be interpreted as a spectral expansion of E{z,p) on the basis of spherical 
harmonics. Near isotropy this expansion converges rapidly, which justifies truncation 
after a few terms. If only the first three terms corresponding to c, m and D are retained, 
a closed system of equations can be derived for these variables by taking moments of the 
conservation equation (2.13) (Baskaran & Marchetti 2009; Saintillan & Shelley 2013). 

In the problem of interest to us here, symmetries dictate that the only non-zero com¬ 
ponents of m and D are mz and Dzz = —2,Dxx = —‘^Dyy in the absence of flow. When 
a flow is applied in the y direction, my and Dyz = Dzy are also expected to become 
non-zero, and Dyy need no longer be equal to D^x- The governing equations for these 
variables can be obtained as 
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(2.19) 

+ (^ + 0”"^ " -^PefSiz)my, (2.20) 

- 2ylPe^+ TOj^ = ‘^PefS{z)m^, (2-21) 

= ^Pe/5(^)P,„ (2.22) 

= -^-PejS{z)Dy^, (2.23) 


Ipe,^ _ 2^Pe2^ + 3P,, = Pe/P(z) (^^c+ . (2.24) 


No equation is needed for D^x, which can simply be deduced from Dyy and D^z using 
the tracelessness of D. In each of these equations, the first term on the left-hand side 
arises due to self-propulsion, the second term captures translational diffusion, and the 
third term rotational diffusion. Terms on the right-hand side arise from the externally 
applied pressure-driven flow and vanish in the absence of flow (Pe/ = 0). 

Boundary conditions for these variables are also readily obtained by taking moments 
of equation (2.14), yielding 


dc 1 
dz 2APes 

dlTlz 1 


niz 


dz 

dDz:z 

d0 


2APes 

2 


Dz 


' 3*^ 
dP,, 


dmj, 

”d7 


1 


15^Pe 


-niz 


dz 


15APe 


2APes^^^' 
dP„ 

-niz 


1 


dz 


lOylPe, 




(2.25) 

(2.26) 
(2.27) 


all to be enforced at z = ±1. For symmetry reasons, we expect c, roy, Dyy, D^z to be even 
functions of z, whereas niz and Dyz are expected to be odd functions. While we consider 
rotational diffusion as the only orientation decorrelation mechanism in this work, all the 
derivations shown here can be easily modified to account for run-and-tumble dynamics 
instead by modifying numerical prefactors in the third terms on the left-hand sides of 
equations (2.22)-(2.24). 

Integrating equation (2.19) and making use of the boundary condition (2.25) easily 
shows that (2.19) can be replaced by 


m. 


- 2APe,— = 0 


dc 

ck 


(2.28) 


at every point in the channel, underlining the direct relation between transverse polariza¬ 
tion and concentration gradients. We also note that the normalization condition (2.12) 
on the distribution function translates into an integral condition on the concentration 
Held expressing conservation of the total particle number: 


c(z) dz = 2. 


(2.29) 


As we discuss next, solution of the system (2.19)-(2.24) subject to the boundary con¬ 
ditions (2.25)-(2.27) and to the integral constraint (2.29) is possible under certain as¬ 
sumptions, and provides results that are in excellent quantitative agreement with the full 
numerical solution of the Smoluchowski equation (2.13) over a wide range of values of 
the Peclet numbers. 
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3. Equilibrium distributions in the absence of flow 

We first analyze the case of no external flow {Pe / = 0), where we expect the boundary 
condition (2.14) to lead to near-wall accumulation and polarization as a result of self¬ 
propulsion. In this case, the full governing equation (2.13) simplifies to 

Pe, (^cos0— - = -Vp, (3.1) 

subject to condition (2.14) at the walls. We note some interesting mathematical properties 
of these equations. First, taking the cross-sectional average of equation (3.1) yields 



which implies that the gap-averaged orientation distribution is isotropic in the absence 
of flow. Using the conservation constraint (2.29), we obtain 

= S' <“> 

which also implies that the first and higher-order moments all average to zero across the 
channel width when there is no flow. 

It is also easily seen that the uniform and isotropic distribution E = 1/An is an 
exact solution of equation (3.1) for all parameter values, though it violates the boundary 
condition (2.14) when A ^ oo. Inspection of the equations shows that, in the limit of 
APcs = dt/2Vs —)■ 0, there is a loss of the higher derivative in both the governing equation 
and the boundary condition. This singular limit suggests the existence of an accumulation 
layer near the channel walls where the distribution departs from the uniform isotropic 
state. Inside this boundary layer, the effects of self-propulsion must be balanced by 
translational diffusion, notwithstanding the small value of APcg ■ Rescaling the governing 
equation inside the boundary layer, however, does not lead to analytical simplifications 
for finite A, so we turn to the simplified moment equations for further characterization of 
particle distributions near the walls in §3.1, where a simple analytical solution is derived 
together with a scaling for the thickness of the accumulation layer. We then describe how 
the limits of strong and weak propulsion can be addressed using asymptotic expansions 
in §3.2 and §3.3. 


3.1. Theory based on moment equations 

In the absence of flow, the moment equations derived in §2.4 only involve c, ruz and Dzz, 
and simplify to: 

dc 

- 2ylPe,—= 0, (3.4) 

az 

Pe.l^-2.1Pe!^ + (T + l)m.=0, (3.5) 

+ (3.6) 

subject to the integral constraint (2.29) and to the boundary conditions 


dmz 

“dT 


1 

2APes 



1 

3^ 


dD 


ZZ 


dz 


15 APes^" 


at z = ±1. 


(3.7) 


Using this set of equations, we first proceed to derive a relation between the values 
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of the concentration and wall-normal polarization at the boundaries. First, we integrate 
equation (3.6) across the channel width and use the second boundary condition in (3.7) 
to arrive at 

J D,,{z)dz = 0. (3.8) 

Now, combining equations (3.4) and (3.5), integrating from z to 1 and making use of the 
first boundary condition gives 

D,,-2APes^ + ^^^c = 2Ac{l). (3.9) 

This relation can be integrated once more across the channel width. Using condition (3.8) 
together with the parity properties of c and m^, this simplifies to 

ci±l)=(^l+-^^TPesm,{±l), (3.10) 


providing a simple relation between concentration and polarization at the walls. Inserting 
this relation into the first condition in (3.7) yields a new set of boundary conditions that 
does not involve the concentration: 


dm^ 1 / 6yl-l-l\ 1 

dz 2APes \ 18yl ) ^ 6A dz 


15APes 


at z = ±l. (3.11) 


Equations (3.5)-(3.6), together with these boundary conditions, form a coupled system 
of second-order linear ordinary differential equations for niz and D^z that can be solved 
analytically. Once these variables are known, the concentration profile is easily obtained 
from the polarization by integration of (3.4) along with condition (3.10). 

Solving these equations yields complicated expressions for c, and Dzz that are 
omitted here for brevity. The prohles, which are illustrated in hgure 2 and will be dis¬ 
cussed in more detail below, reveal one important finding: while a significant wall-normal 
polarization exists in the near-wall region, nematic alignment is relatively weak through¬ 
out the channel for A > 0.1. This suggests seeking a yet simpler solution that neglects 
nematic order altogether. If the moment expansion (2.18) is truncated after two terms, 
the equations for c and rUz simplify to 


niz — 2APes 


dc 

(U 


= 0 , 


-2APel 


dfmz 

dz^ 


-b 




rUz 


= 0 , 


subject to the conditions 


druz 

“d7 


c 

QAPcs 


at z = ±1, 


and 


J c(z) dz = 2. 


(3.12) 


(3.13) 


Solving these equations is straightforward and provides elegant expressions for the con¬ 
centration and polarization profiles: 


where 


B [6A cosh B + cosh Bz] 

6AB cosh B + sinh B ’ 

(3.14) 

GAPCsB^ sinhBz 

TfL ( z) — 

^ 3 (GtIB coshB-b sinhi?) ’ 

(3.15) 


(3.16) 
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defines the dimensionless decay length of the excess concentration at the walls. In dimen¬ 
sional terms, this decay length is given by B~^H = f?a\/3/(l + &A) where £a = dtjVa- 
In the limit of strong propulsion {A <C 1), it simplifies to \fZia- In the limit of weak 
propulsion {A 1), it becomes where id = \fdijdr is a purely diffusive length 

scale. For Brownian particles, id is typically of the order of the particle size L, though 
this may not be the case for active particles subject to athermal sources of noise. Next, 
we focus more precisely on these two limits by rescaling the governing equations with 
the appropriate scales identified here. 


3.2. Strong propulsion limit: —)■ 0 

In the limit of small A, the above discussion suggests rescaling the Smoluchowski equation 
using the accumulation length scale ia, yielding 


dE , 2 


subject to the boundary condition 

dE 


dz 


= cos 6E at z = ±H*. 


(3.17) 


(3.18) 


Here, H* = {2APes) ^ is the channel half-height rescaled by the accumulation length 
scale ia- The gap-averaged isotropy constraint is now expressed as 




H* 


(3.19) 


/ Edz = 

J-H- 27T 

The leading-order solution corresponding to yl = 0, which was previously obtained by 
Elgeti & Gompper (2013), is written 

H* cos 9 


dTTsinh (H* cos 0) 


(3.20) 


and it is easily seen that it satisfies zero wall-normal flux pointwise throughout the chan¬ 
nel. In particular, it shows that wall accumulation is possible even in the absence of 
rotational diffusion and is simply a result of a coupling between self-propulsion, transla¬ 
tional diffusion and confinement. This solution can then be corrected to order 0{A) by 
solving the first-order inhomogeneous equation 

rld/A) p 

cos6E^^\z,e)-=Vl tf'(°)(z,6»)dz. (3.21) 

OZ J-H- 


subject to boundary condition (3.23). An exact analytical solution to this equation can 
again be obtained but is cumbersome and omitted here for brevity. 


3.3. Weak propulsion limit: A —> oo 


In the limit of large A, the Smoluchowski equation is rescaled using the diffusive length 
scale id as 


subject to 




cos 



dz'^ 



(3.22) 


dE 1 

cos 0E at z = ±H\ (3.23) 

dz ^ 


where = {2y/APes) The leading-order solution in the limit of A —)■ oo is uniform 
and isotropic and corresponds to the case of a passive particle: = l/An. It can 
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Figure 2. (Color online) Equilibrium distributions in the absence of flow and for various swim¬ 
ming Peclet numbers Pcs (with A = 1/6), obtained by numerical solution of equation (2.13) 
using finite volumes: (a) concentration c, (6) wall-normal polarization m^, and (c) wall-normal 
nematic order parameter 


be corrected asymptotically using a regular perturbation expansion in powers of A 

T {z, 0) = {z, 0) + yl-i /2 ^( 1 ) ^-1 ^( 2 ) _ _ (3 24 ) 

Recursively solving for higher-order terms yields 


(z,6») 

(z,6») 


3 

-—--—— sinh(v^ 2 :) cos 0, 

AnV2cosh{V2W) 

1 cosh-\/2z tanh(-\/2i7i) cosh(-\/6z) / 1\ 

— ---- I COS 0 — TT I i 

15 cosh(-\/2i/t) 5-\/3 sinh(-\/6iJt) V 3/ 


(3.25) 

(3.26) 


which both satisfy the appropriate boundary conditions. Quite remarkably, it can be 
seen that successive terms in the expansion (3.24) correspond to successive orientational 
moments of the distribution function in equation (2.18), with and describing 
the polarization and nematic order, respectively. 


3.4. Numerical results and discussion 

Figure 2 shows the full numerical solution for the concentration c, wall-normal polar¬ 
ization mz and nematic order parameter Dzz obtained by finite-volume solution of the 
Smoluchowski equation (2.13) as described in Appendix B. Here, we fix the value of A 
and focus on the effect of Pcs, which an inverse measure of confinement. The concentra¬ 
tion profiles shown in figure 2(a) exhibit significant accumulation of particles near the 
boundaries, especially at low values of Peg. As anticipated, this accumulation is accom¬ 
panied by polarization towards the boundaries as a direct consequence of the boundary 
condition (2.25), as well as by a weak nematic alignment. As Pcg increases, the spa¬ 
tial heterogeneity and anisotropy near the walls progressively extend through the entire 
channel as the two boundary layers thicken and eventually merge. Further increase in 
the swimming Peclet number leads to a flattening of the profiles, which is especially 
significant when Peg > 1. This flattening is a direct consequence of the scaling of trans¬ 
lational diffusion with Pe^ in equation (2.13), causing it to overwhelm self-propulsion 
which scales with Peg. The influence of A is illustrated in figure 3, where it is seen to 
be similar to that of Peg: increasing A leads to a thickening of the boundary layers and 
flattening of the concentration profiles, again due to the scaling of translational diffusion 
with A in equation (2.13). 

The finite-volume numerical solution of the full conservation equation (2.13) is in 
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Figure 3. (Color online) Equilibrium distributions in the absence of flow and for various values 
of A (with Pcs = 0.25), obtained by numerical solution of equation (2.13) using finite volumes: 
(a) concentration c, (b) wall-normal polarization niz, and (c) wall-normal nematic order param¬ 
eter Dzz- Solutions based on moment equations are nearly identical, as illustrated in figure 4. 



Pes 


Figure 4. (Color online) The relative rms error for the concentration between the finite-volume 
solution and the two-moment analytical solution (3.14) for different values of A. Solutions based 
on moment equations are nearly identical to the finite-volume solution for sufficiently large 
values of A. 


excellent quantitative agreement with the two- and three-moment approximations derived 
previously, which are not shown in figure 2 as they are nearly indistinguishable over the 
entire channel width as long as > 0.1. The rms error between the two-moment solution 
of equation (3.4) and the finite-volume solution is indeed plotted in hgure 4, where it 
remains below 10“^ for all values of Pcg considered here when A > 0.1. This finding 
may seem quite surprising considering the strong approximation made when truncating 
expansion (2.18) after only two terms, and strongly validates the use of approximate 
moment equations such as (2.19)-(2.24) when modeling active suspensions, at least in 
the absence of flow. For very small values of A, however, nematic alignment at the walls 
becomes significant as seen in figure 3(c), so that the nematic tensor can no longer be 
neglected and the two-moment solution loses its accuracy; in this case, the alternate 
expressions derived in the small A limit in §3.2 can be used instead. 

The influence of PCg on wall accumulation is analyzed more quantitatively in figure 5, 
showing the values of the wall concentration c(±l), the boundary layer thickness 6 defined 
as the distance from the wall where c(l — <5) = 1, and the fraction 6* of particles inside 
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Pe, 


Pes 


Pes 


Figure 5. (Color online) Wall accumulation in the absence of flow as a function of Pes (at 
A — 1/6): (a) concentration c(±l) at the walls; (b) boundary layer thickness 5, defined as the 
distance from the wall where c(l — 5) = 1; (c) fraction S* of particles inside the boundary 
layer, defined as the integral of c{z) over the boundary layer thickness. The solid line shows the 
theoretical prediction based on the two-moment solution (3.14), and symbols show full numerical 
results using finite volumes. 


the boundary layer defined as 

S* = f c{z) dz. 
Jl-5 


(3.27) 


Analytical expressions for these quantities can be derived from the two-moment solution 
(3.14). In particular, the boundary layer thickness is obtained as 


6{Pes) = 1 - ^log 


sinhi? / sinhi? 
B B 



(3.28) 


which has the two limits 


lim 6(Pes) = 0 and lim 6{Pes) = 1- 

Pbs^O Pbs^oo ^3 

Similarly, the fraction of particles inside the boundary layer is given by 

6AB (1 — 6) coshi? + sinh [B (1 — 5)] 


S*iPes) = l- 


6AB cosh B + sinh B 


(3.29) 


(3.30) 


and has the same limits as 5{Pes) when Pcg —t 0 and oo. 

As shown in figure 5(a), the wall concentration reaches its maximum in the limit of 
Pcs —>■ 0, and steadily decreases towards 1 as Pcs increases due to the smoothing effect of 
translational diffusion. This is accompanied by an increase in the boundary layer thickness 
6, which asymptotes at high values of Pcg- The fraction 5* of particles near the walls 
shows a similar trend, but interestingly also exhibits a weak maximum for Pe^ ~ 1.135 
when wall accumulation due to self-propulsion and translational diffusion are of similar 
magnitudes; at this value of Peg, 6* « 0.46 corresponding to nearly half the particles 
being trapped near the walls. As previously observed in figure 4, excellent agreement is 
obtained between the two-moment approximation and the numerical solution of the full 
governing equations. 
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4. Equilibrium distributions and transport in flow 

4.1. Weak-flow limit: regular asymptotic expansion 

We now proceed to analyze the effects of an external pressure-driven flow, first focusing 
on the case of a weak flow for which Pep <C 1. Since the parameter A is fixed for a given 
type of swimmers, we keep it constant in the rest of the paper and focus on the effects 
of Pcs and Pej. The form of the governing equations suggests seeking an approximate 
solution as a regular expansion of the moments of the distribution function in powers 
of Pcf. The leading-order O(Pe^) solution corresponding to the absence of flow was 
previously calculated in §3. It is henceforth denoted by and we recall 

that = Dy^J = 0. Inspection of the moment equations (2.I9)-(2.24) reveals that the 
interaction of the applied shear profile S{z) with this leading-order solution perturbs my 
and Dyz at order O(Pef). On the other hand, c, m^, D^z and Dyy are only perturbed 
by the flow at order 0{Pe^j) due to its interaction with my and Dyz. Based on these 


observations, we expand the solution as 

c{z) = c^°\z) + Pe) c(2)(z) -P 0{Pe)), (4.1) 

mz{z) = m'^^\z) + Pe^fm^^\z) P 0{Pe^f), (4.2) 

Dzz{z) = DfXz) + PejD^^Xz) + 0{Pej), (4.3) 

Dyyiz) = Df^{z) + Pe}D^;fJ{z) + OiPej), (4.4) 

my{z) = Pefm^^\z) + 0{Pe)), (4.5) 

Dyz{z) = PefDg^z) + 0{Pe}). (4.6) 


We focus here on determining the leading-order corrections to my and Dyz, which capture 
streamwise polarization and nematic alignment with the applied shear, respectively. The 
0{Pef) moment equations are written 


Pe 


dD 


( 1 ) 


dz 


- 2APe 


,d2 


m 


( 1 ) 




■ ml}^ = --S{z)m 
^ 5 


( 0 ) 


Pe., dm 


( 1 ) 


dz 


- 2ylPe: 


dz^ 


+ 3PW = S{z) , 


subject to boundary conditions 

dmy ^ _ 1 dDy) 

dz ~ 2APes ’ dz 


_^ 

loffPe, y 


at z = ±1. 


(4.7) 

(4.8) 


(4.9) 


Note that the forcing terms on the right-hand sides of equations (4.7)-(4.8) are known 
and capture the interaction of the local shear rate S'(z) with the equilibrium distributions 
in the absence of flow. 

A numerical solution of equations (4.7)-(4.9) is plotted in figure 6 for different values 
of Peg. At low values of the swimming Peclet number, figure 6(a) shows an upstream 
polarization (ruy < 0) near the boundaries, and a downstream polarization (ruy > 0) near 
the center of the channel. The upstream polarization, which has previously been observed 
in both experiments and simulations and is at the origin of the well-known phenomenon 
of upstream swimming, is a simple and direct consequence of the shear rotation of the 
particles near the wall, which tend to point towards the walls in the absence of flow as 
explained in §3. This interaction is encapsulated in the right-hand side in equation (4.7). 
The downstream polarization near the centerline is a more subtle effect arising from self¬ 
propulsion through the first term on the left-hand side of (4.7). As Pe^ increases and the 
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2 2 


Figure 6. (Color online) Effect of a weak applied flow: leading-order 0{Pef) corrections of 
(a) streamwise polarization rriy and (b) shear nematic alignment Dy^ for different values of the 
swimming Peclet number, obtained by numerical solution of equations (4.7)-(4.9). 


boundary layers thicken, upstream swimming becomes weaker near the boundaries due 
to the weaker wall-normal polarization there; however, my is also observed to become 
negative across the entire channel due to the thickening of the polarized boundary layers 
into the bulk of the channel as previously shown in figure 2(6). 

The mean streamwise swimming velocity Vy of the active particles with respect to the 
imposed flow can be defined in terms of the polarization as 

Vy = ^J Pes my{z) dz = J (z) dz = PcsPefTn^y^ (4.10) 


An expression for can be derived based on the moment equations. We first take the 
cross-sectional average of equation (4.7) and use the first boundary condition to obtain 


m 


(1) = -- 
V 5 



zm<f\z) dz. 


(4.11) 


Since is an odd function of z with m^^\z) > 0 for z > 0, the integrand on the 
right-hand side is always positive across the channel, and therefore the mean upstream 
polarization is negative: < 0. This also implies that Vy < 0, i.e., there is a net 

upstream flux of particles against the mean flow for all values of A and Peg in the weak 
flow limit. Using equation (3.5) for m^^\z), we can rewrite the right-hand side as 


m^y^ = — 


+ 1) 


2 f d'^ra 


(0) 


dz — Peg 


dD^J 

dz 


dz 


(4.12) 


After integration by parts and application of the boundary condition on (z) together 
with equation (3.8), this simplifies to 


m(l) = — 


2Peg 




(.( 0 ) ( 1 ) _ QAPcg ( 1 ) 


(4.13) 


Recalling that c^°)(l) and mi°^(l) are related via equation (3.10), we obtain two expres¬ 
sions for the mean streamwise swimming velocity in terms of either the concentration or 
wall-normal polarization at the top wall in the absence of flow: 


— AA ^ 

Vy = -—PelPef 


;(o)(i)_i pg2pg^ 1 - dAPegm^°\l) 

J ID L 


(4.14) 
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Since the concentration at the wall in the absence of flow always exceeds the mean when 
PCs > 0, equation (4.14) again confirms that Vy < O.If we further make use of the 
simplified two-moment analytical solution (3.14) for the concentration profile, we arrive 
at a simple expression for the mean upstream velocity in terms of the swimming and flow 
Peclet numbers: 


V 


y ~ 


4.1 

T 


PelPef 


B cosh B — sinh B 
QAB cosh B + sinh B 


(4.15) 


This simple analytical prediction for Vy will be tested against numerical simulations at 
arbitrary Ptf in §4.2, where it will be shown to provide an excellent estimate for the 
swimming flux up to Ptf « 2. 

The effects of the external flow on nematic alignment are also illustrated in figure &{d), 
where Dyz is found to vary almost linearly across the channel width and has the same sign 
as the external shear rate profile S{z). The right-hand side in equation (4.8) provides a 
simple explanation for these findings, where we see that shear nematic alignment results 
primarily from the interaction of the flow with the concentration profile and with the 
wall-normal nematic alignment. As PCg increases, shear nematic alignment decreases due 
to the decrease in c and D^z inside the boundary layers as seen in figures 2(a) and (c), 
and to self-propulsion through the first term on the left-hand side of equation (4.8). 


4.2. Strong-flow limit: scaling analysis 

As we shall see in §4.3 and figure 7, the regime of high flow Peclet number is also quite 
interesting as it can result in a depletion near the channel centerline surrounded by regions 
where particles become trapped. The thickness of this depletion region will be found to 
decrease with increasing flow strength, suggesting the presence of another boundary layer 
near 2 : = 0 in the limit of Peg ^ 1. Insight into this regime can be gained by analyzing 
the behavior of the governing equation (2.13) for Pcf 1 and Pcg <C 1. If the swimming 
Peclet number is low, the wall boundary layers are very thin and have negligible impact 
on the dynamics in the bulk of the channel. Inspection of equation (2.13) suggests that, 
in the outer region away from both the channel walls and the centerline, the dominant 
balance is between shear alignment and rotational diffusion: 

^ S{z)Vp • [cos9{l -pp) • yif] « 


In this region, the concentration is expected to be nearly uniform, and the particle orien¬ 
tation distribution is primarily nematic as a result of the competition between the local 
shear rate and rotational diffusion (as would occur in a passive rod suspension). This 
corresponds to the shear-trapping region where cross-streamline migration is very weak 
due to the strong alignment with the flow. 

However, as we move closer and closer to the centerline, the local shear rate decreases, 
causing a concomitant decrease in shear alignment and increase in cross-streamline mi¬ 
gration due to self-propulsion. This transition corresponds to the edge of the central 
boundary layer from which particles are depleted, and the position of this transition 
region (or half-thickness of the depletion layer) can be estimated by balancing the mag¬ 
nitudes of the terms describing self-propulsion and shear alignment in equation (2.13): 


Pes 

Sd 



(4.17) 


from which we find 


Sd - 


(4.18) 
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where the prefactor C is a numerical constant and where we have defined 

^ Pe, ^ Fs 
^ Pef 27 „iJ' 


(4.19) 


The dimensionless group x can be interpreted as the ratio of the timescale 7 “^ it takes a 
particle to align with the flow over the characteristic timescale 2H/Vs it takes it to swim 
across the channel width: if x is small, particles align with the flow much faster than 
they can cross the channel, leading to significant shear-trapping; on the other hand, if x 
is large, particles cross the channel much faster than they align with the flow and shear¬ 
trapping does not occur. As we show in Appendix C, this scaling for 6d can indeed also 
be derived by considering the individual trajectories of deterministic swimmers released 
from the centerline, which can be shown to become trapped at a distance of the order of 
(5£). It will also be shown to agree quite well with numerical results in §4.3, where we will 
find that Sd ~ 2.404y'x provides an excellent estimate for the thickness of the depletion 
layer when Pcs < 0.25 and Pef > 50. 

To gain further understanding of the effect of shear rate on the intensity of depletion, 
we rescale lengths by (5 d inside the central boundary layer to rewrite the governing 
equation (2.13) as 


P dT 

— cos 6 »-- 2 A 

C oz 


7-2 q2^ 


CP 


-zV 


p ■ ^cose{l -pp) ■yT]= 


(4.20) 


where the dimensionless group P = y^PCsPcf emerges as the most significant pa¬ 
rameter governing the profile of the depletion layer. Unsurprisingly, we find that self¬ 
propulsion and shear rotation have the same magnitude upon rescaling. In this region, 
self-propulsion, which scales with P, has the effect of enhancing depletion by driving par¬ 
ticles away from the centerline; this competes against translational diffusion, scaling with 
P^, which has the effect of smoothing concentration gradients and thus hampers deple¬ 
tion. This suggests the following dependence of the concentration profile on Pcf. As flow 
strength is increased from small values, the depletion layer forms and continually narrows 
according to equation (4.18) for Sd- As long as T < 1 , self-propulsion dominates trans¬ 
lational diffusion and increasing Pcf (and therefore T) enhances depletion. This trend 
reverses when P ~ 0 ( 1 ), when translational diffusion starts to overcome self-propulsion, 
leading to a subsequent decrease in the strength of depletion for T > 1. This qualitative 
explanation for the non-monotonic dependence of the strength of depletion upon P (and 
hence upon the mean shear rate of the imposed Poiseuille flow) is consistent with the 
experimental observations of Rusconi et al. (2014), and is also borne out by numerical 
solutions of the governing equations as we describe next. 


4.3. Arbitrary flow strengths: finite-volume calculations and discussion 

We now test and extend the key predictions from the weak-flow asymptotics and strong- 
flow scaling analysis from the preceding sections by performing finite-volume numerical 
simulations of the governing equation (2.13) for arbitrary values of Pcs and Pef using 
the algorithm of Appendix C. Typical concentration profiles are illustrated in figure 7 for 
various values of Pe/, and for the two values of Pcs = 0.25 and 1.0 corresponding to cases 
where wall accumulation in the absence of flow is strong and weak, respectively. In both 
cases, the leading effect of the external flow on c is to decrease wall accumulation. This 
trend is easily understood as a result of the alignment of the particles with the flow, which 
reduces wall-normal polarization and thereby hinders accumulation. This decrease in 
accumulation also results in a net increase in the concentration in the central parts of the 
channel and in the flattening of the profiles in the strong-flow limit. When Peg is small as 
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Figure 7. (Color online) Equilibrium concentration profiles (at A = 1/6) for (a) Pea = 0.25 
(strong wall accumulation) and (b) Pcs = 1.0 (weak accumulation) and for various values of the 
flow Peclet number Pej, obtained by finite-volume solution of the governing equation (2.13). 

in figure 7(a), a depletion layer is also observed to form near the channel centerline and to 
progressively narrow with increasing Pe /, in agreement with the theoretical predictions of 
§4.2. At high values of Pe/, the three distinct regions identified in §4.2 (wall accumulation, 
shear-trapping, and centerline depletion) in fact become clearly visible. However, if the 
swimming Peclet number is increased to Peg = 1.0 as in figure 7(&), the thickening of the 
wall boundary layers suppresses shear-trapping and depletion at the centerline, leading 
to a nearly uniform concentration profile in the strong flow limit. 

Corresponding profiles for the wall-normal and streamwise polarization are also shown 
in figure 8. As expected, rotation of the particles by the flow causes a decrease in the 
wall-normal polarization, and also results in a non-zero streamwise polarization my as 
previously discussed in §4.1. This streamwise polarization is especially strong in the 
near-wall region where my is negative, indicating upstream swimming. It is significantly 
weaker near the center of the channel, where it is found to be positive for Peg = 0.25 
but remains negative across the entire channel when Peg = 1.0 due to the overlap of the 
two wall boundary layers. 

These trends are made more quantitative in figure 9, showing the dependence of c(±l), 
my(±l) and my(0) on the swimming and flow Peclet numbers. As previously discussed, 
the wall concentration is seen to decrease with increasing flow strength irrespective of the 
value of Pcg, and asymptotically tends to 1 in the strong-flow limit as the concentration 
profiles flatten. Figure 9(6) shows that the streamwise polarization at the walls is always 
negative, which implies that the active particles always swim upstream near the bound¬ 
aries. Interestingly, we find that there is maximum upstream swimming at Pe/ « 10, and 
the upstream motion is reduced at higher values of the flow Peclet number. The stream- 
wise polarization at the channel centerline shows complex trends as shown in figure 9(c). 
As predicted by the weak-flow asymptotic analysis of §4.1, TOy(O) is found to be positive 
for low values of Pcg and negative for high values of Pcg. Its absolute value increases 
with flow strength in both cases up to Pe/ « 10, beyond which further increasing flow 
strength reduces the polarization. The decrease in both mj,(±l) and my(0} at high Pe/ 
is a likely consequence of the dominant effect of the shear alignment term in equation 
(2.13), which promotes nematic rather than polar order. 

The dependence of the average streamwise swimming velocity Vy defined in equa¬ 
tion (4.10) on both Peclet numbers is shown in figure 10, where numerical results are 
compared to the weak-flow theoretical prediction of equation (4.14). Consistent with fig¬ 
ure 9(6) for the streamwise polarization at the walls, we find that Vy < 0, and that \Vy\ 
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Figure 8. (Color online) Equilibrium streamwise and wall-normal polarization profiles (at 
A = 1/6) for (a)-(c) Pca — 0.25 and (6)-(d) Pea = 1.0 and for various values of the flow 
Peclet number Pcf, obtained by finite-volume solution of the governing equation (2.13). The 
streamwise polarization niy is shown on the top row (a)-(fc), and the wall-normal polarization 
rriz on the bottom row (c)-(d). 



Pef Pef Pef 


Figure 9. (Color online) Effect of swimming and flow Peclet numbers on: (a) wall concentration 
c(±l), (6) streamwise polarization my{±l) at the channel walls, and (c) streamwise polarization 
my{0) at the channel centerline. 


first increases nearly linearly with Pef in agreement with the predictions of §4.1. This 
increase persists up to Pef « 10, beyond which \Vy\ starts decreasing again. Excellent 
quantitative agreement is found with equation (4.14) for Pcf ^ 2.0. This is confirmed in 
hgure 10(6), showing the dependence of \Vy\/Pef on swimming Peclet number: the up¬ 
stream velocity is found to increase with Pcg, primarily as a result of the corresponding 
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Figure 10. (Color online) (a) Magnitude of the average upstream swimming velocity \Vy\ as 
a function of Pcf for different values of Pcs (at A = 1/6), and (6) dependence of \Vy\/Pef on 
Pes for different values of Pcf. Symbols show finite-volume numerical simulations, and dotted 
lines show the theoretical prediction of equation (4.14). 


increase in swimming speed of individual particles, and a collapse of all the curves onto 
the theoretical prediction of equation (4.14) is observed when Pef < 2.0. 

As seen in figure 7(a), shear-trapping and centerline depletion are observed in the 
central portion of the channel at high flow Peclet number if Peg is sufficiently low. This 
is illustrated more clearly in figure 11, where concentration and wall-normal polarization 
profiles are shown in the central portion of the channel for various values of the flow Peclet 
number and for Pcg = 0.125. This value was chosen to match the experiments of Rusconi 
et al. (2014), where the following parameters were reported: Vg = 50 /im, = 1 s“^, and 
2H = 400 /im. As seen in figure 11(a), increasing Pe/ from zero first results in a decrease 
in the concentration at the centerline, corresponding to the formation of the depletion 
layer. As the concentration decreases, the width of the depletion layer is also found to 
decrease. This trend continues up to Pef « 20, above which the concentration at the 
centerline starts increasing again, even though the depletion layer keeps narrowing. These 
trends are in very good agreement with the experiments of Rusconi et al. (2014), who 
also reported a non-monotonic dependence of the strength of depletion on shear rate; 
in fact, the profiles shown in figure (11) are very similar to the experimental profiles at 
equivalent values of Pef. The trends on the concentration profile are easily understood 
based on figure ll(fe) for the wall-normal polarization, which reflects the net swimming 
velocity across the channel and provides insight into cross-streamline migration. Indeed, 
the polarization profiles exhibit peaks on both sides of the depletion layer, corresponding 
to a strong migration away from the center. These peaks increase in magnitude and also 
shift towards the centerline as flow strength increases and the depletion layer narrows. 
Beyond those peaks, quickly decays to zero where the concentration profiles plateau 
in accordance with equation (2.28) and shear-trapping of the particles takes place. 

These trends are tested more quantitatively against the strong-flow scaling analysis of 
§4.2 in figure 12. We first define the thickness of the depletion layer as the distance 
from the centerline where reaches its maximum, when such a maximum exists. Based 
on the analysis of §4.2, we expect 5d to scale linearly with ^/x = Peg/Pef in strong 
flows, and this is indeed confirmed in figure 12(a). We find that 5d can only be defined 
when 0.16 or Pef > 40Pes, which corresponds to the shear-trapping regime. Best 

agreement with the scaling prediction is obtained in the low Peg and high Pef limit, and a 
linear least-square fit to the data for Peg < 0.25 and Pef > 50 shows that 5d « 2.404^/%. 
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^ 2: 


Figure 11. (Color online) (o) Concentration profiles in the central portion of the channel for 
Pcs = 0.125 and various values of the flow Peclet number Pcf, obtained by finite-volume 
solution of equation (2.13). (6) Corresponding profiles of the wall-normal polarization m^. 



Figure 12. (Color online) (a) Depletion layer thickness So, defined as the distance from 
the centerline where the wall-normal polarization reaches its maximum, as a function of 
.^/X = \JPcsjPcf. (b) Depletion index Ad defined in equation (4.21) as a function of 
P = y/PCsPCf. 


As PCs increases, the numerical results depart from this prediction, primarily due to the 
thickening of the wall boundary layers which causes them to interact with the parts of the 
channel where shear-trapping and depletion occur. We further quantify the shape of the 
depletion layer by introducing a depletion index Ad measuring the amount of particles 
depleted from the center due to trapping in high-shear regions: 

Ad = [ c{z)dz - Sdc{6d)- (4-21) 

Jo 

As we argued in §4.2 based on equation (4.20), the shape of the depletion layer is expected 
to depend upon P = y/PegPef, and indeed the numerical data for the depletion index 
for various values of Pcs and Pcf is found to collapse onto a master curve when plotted 
vs P in figure 12(6). In agreement with the trends observed in figure 11(a), the depletion 
index shows a non-monotonic dependence on A, with maximum depletion occurring for 

r « 2. 

The dynamics in the limits of Pcs 1 and Pe / 1 are summarized schematically in 

figure 13, where the channel can be roughly divided into three distinct regions. Region 
(A), with thickness 5 ^ APcg, abuts the channel wall and is characterized by wall accu- 
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Figure 13. Schematic summary of the dynamics in the limits of Pcs ^ 1 and Pe/ 1. The 
channel can be roughly divided into three regions: (A) near the walls, particles accumulate in a 
boundary layer of thickness <5 ~ APcs', (B) away from the walls and centerline, strong nematic 
alignment by the flow leads to shear-trapping and a nearly uniform concentration profile; (C) 
near the centerline, particle propulsion leads to a depletion layer of thickness 5d ~ P- The 
diagram only shows the left half of the channel 2 G [—1,0]; the corresponding diagram in the 
other half can be obtained by symmetry and by noting that is an even function of z, whereas 
TUy and Dy^ are both odd functions. 


mulation and a net polarization towards the wall. These effects occur even in the absence 
of flow, and are in fact mitigated by the flow which tends to decrease the wall concentra¬ 
tion and rotate particles to induce upstream polarization. Away from both the wall and 
the channel centerline is region (B), where the concentration profile is nearly uniform and 
shear trapping occurs: here, polarization is weak but there is a strong nematic alignment 
of the particles due to the applied shear. The local shear rate decreases in magnitude as 
we approach the centerline and enter region (C), which has a characteristic thickness of 
5d ~ y^PCs/Pef. in this region, particles are depleted due to a net polarization towards 
the walls, which drives migration away from the center but is counterbalanced by trans¬ 
lational diffusion. Increasing Pcg causes both regions (A) and (C) to widen, up to a point 
where they merge and the three regions can no longer be distinguished. Increasing Pcf, 
on the other hand, tends to weaken wall accumulation but does not change the thickness 
of region (A), while it also causes the narrowing of region (C). 


5. Discussion 


5.1. Summary of main results 

We have used a combination of theory and numerical simulations to analyze the dis¬ 
tributions and transport properties of an infinitely dilute suspension of self-propelled 
particles confined between two parallel flat plates, both in quiescent conditions and un¬ 
der an imposed pressure-driven flow. Our analysis focused on incorporating the effects of 
confinement within the kinetic theory framework previously developed by Saintillan & 
Shelley (2008a), which is based on a Smoluchowski equation for the distribution of the 
active particle positions and orientations. In particular, we demonstrated that prescrib¬ 
ing a zero-normal-flux condition on the particle distribution function at the boundaries 
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captures several key features reported in experiments on dilute active suspensions un¬ 
der confinement. We presented a finite-volume algorithm for the numerical solution of 
the Smoluchowski equation, which allows for an easy implementation of the boundary 
conditions, and also developed a simpler system of equations for the orientational mo¬ 
ments of the distribution function, which enabled us to perform analytical calculations 
in the absence of flow and under a weak imposed flow. An asymptotic scaling analysis 
was also performed on the full Smoluchowski equation under strong flow. The numerical 
simulation data was used to test and further understand the analytical calculations and 
predictions. 

We hrst considered the dynamics in the absence of flow. In this case, the governing 
equations involve a swimming Peclet number Pe^, which is the ratio of the persistence 
length of swimmer trajectories to the channel height, as well as a parameter A that is 
hxed for a given swimmer type and whose inverse measures the strength of propulsion. 
In the limit of wide channels, the channel can be divided into two regions: a near-wall 
accumulation region where the particles tend to concentrate and have a net polarization 
towards the wall, and a bulk region away from the walls where the distribution is nearly 
uniform and isotropic. Asymptotic expressions for the full distribution function were also 
derived as series in powers of A in the weak and strong propulsion limits. In particular, it 
was shown that the characteristic thickness of the accumulation layer scales with dt/Vs in 
the strong propulsion limit (A <C 1), and with \fdt]~dr in the weak propulsion limit {A ^ 
1). For finite values of A, analytical expressions for the concentration and polarization 
prohles were obtained by solving the moment equations and displayed excellent agreement 
with the finite-volume numerical simulation of the full distribution function for a wide 
range of values of the swimming Peclet number so long as A > 0.1. Based on these 
results, we proposed and validated a simple mechanism for wall accumulation, where 
the presence of the wall breaks the polar symmetry of the active particles and leads 
to sorting of orientations. This mechanism differs from previous explanations based on 
hydrodynamic interactions or surface alignment due to collisions, and led us to conclude 
that both pusher and puller particle suspensions will exhibit similar wall accumulation in 
the dilute limit. Hydrodynamic and surface alignment interactions are, however, expected 
to quantitatively affect the profiles in more concentrated systems and to lead to different 
distributions for pusher and puller particles. 

Next, we analyzed the effects of an imposed pressure-driven flow. When a flow is 
applied on the suspension, the physics is now governed by three dimensionless groups: 
the swimming Peclet number Pcs and parameter A introduced above, as well as a flow 
Peclet number Pcf comparing the imposed shear rate to rotational diffusion. In the weak 
flow limit, we calculated the leading-order corrections of the streamwise polarization and 
shear nematic alignment due to the flow and showed that near-wall upstream swimming 
is a consequence of shear rotation of the particles inside the accumulation layer near the 
walls. We derived an analytical expression for the average upstream swimming velocity of 
the active particles relative to the imposed flow, which was compared against numerical 
simulations and provides an excellent estimate for Pe/ < 2. In the strong flow limit, we 
developed a scaling analysis to show that when Peg ^ 1 and Pe / 1 the channel can 

be roughly divided into three regions: the near-wall accumulation region with thickness 
S ~ APcs, a depletion region near the centerline with thickness 5d ^ P = \/Pgs/P^ f, 
and a shear-trapping region away from the wall and centerline where the concentration 
is nearly uniform and particle alignment is primarily nematic. The extent of the central 
depletion shows a non-monotonic variation with flow strength, with a maximum depletion 
occurring at a critical flow strength such that P ~ 0(1). 
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5.2. Discussion and comparison to previous works 


The phenomena analyzed in this study have received considerable attention in experi¬ 
ments as well as other models and simulations, so we compare and contrast them here to 
these prior works. As mentioned in the introduction, the wall accumulation predicted by 
our model in the absence of flow is well known in experiments on bacterial suspensions, 
where accumulation layers of « 1 to 50 ^m are typically reported (Berke et al. 2008; 
Li & Tang 2009; Li et al. 2011; Gachelin et al. 2014), with increases in concentration 
of up to 50 times the bulk density very close to the wall (Li et al. 2011). Such high 
concentrations at the walls are consistent with our numerical results of figure 3, which 
predict high values of c(±l) in the strong-propulsion limit of A <C 1 relevant to bacteria. 
Indeed, a rough estimate for E. coli provides A « 0.01, though it is difficult to precisely 
measure dt in experiments since long-time mean-square displacements are dominated by 
Taylor dispersion. This strong accumulation is also consistently observed in simulations 
(Hernandez-Ortiz et al. 2005; Nash et al. 2010; Costanzo et al. 2012; Elgeti & Gompper 
2013; Lushi et al. 2014; Li & Ardekani 2014), which also exhibit the preferential align¬ 
ment of the swimmers towards the wall that our model predicts. A similar alignment 
has also been reported in a few experiments (Drescher et al. 2011; Lushi et al. 2014), 
though detailed observations of swimming micro-organisms near walls has also revealed 
complex complex scattering dynamics due to the interactions of the flagellar appendages 
with the boundaries (Denissenko et al. 2012; Kantsler et al. 2013). These observations 
seem to contradict mechanisms purely based on Stokes-dipole hydrodynamic interactions 
with the no-slip walls, as these predict reorientation of the cells parallel to the walls in 
the case of pushers (Berke et al. 2008). Rather, they appear to support the prediction 
that accumulation layers derive predominantly from a polarity-sorting mechanism across 
the channel together with a balance of self-propulsion and diffusion at the walls. We note 
that this mechanism was also proposed in the work of Elgeti & Gompper (2013), who 
performed simulations of self-propelled Brownian spheres between two flat plates. Their 
numerical results support the trends described in §3.4 on the effect of confinement as cap¬ 
tured by Peg. Elgeti & Gompper (2013) also wrote down a continuum model that shares 
similarities with ours, which they used to analyze the strong propulsion and narrow gap 
limits. Their conclusions are in agreement with the discussion of §3.2 and §3.3. 

The distributions and dynamics predicted by our theory under imposed flow also agree 
with the bulk of prior studies, both experimental and numerical. The reorientation of 
near-wall swimmers against the flow leading to upstream swimming has been reported 
ubiquitously in many experiments (Hill et al. 2007; Kaya & Koser 2009, 2012; Kantsler 
et al. 2014) and simulations (Nash et al. 2010; Costanzo et al. 2012; Chilukuri et al. 2014), 
with several of these studies proposing similar mechanisms as that described herein, 
namely the shear rotation of the polarized cells near the walls. Quite remarkably, the 
peak in the upstream swimming flux at a critical flow strength visible in the simulation 
data of figure 10(a) was also reported in the experiments of Kantsler et al. (2014). 

The dynamics in strong flows in the central part of the channel has only received 
little attention in previous studies. Our interest in this problem was sparked by the 
recent microfluidic experiments of Rusconi et al. (2014), which were the first to predict 
centerline depletion and shear trapping. Our scaling analysis and numerical results of §4.2 
and §4.3 are in excellent agreement with their observations. In particular, the shape of the 
concentration profiles near the channel centerline obtained in figure 11 are quite similar 
to those shown in figure 2(a) of their paper. Further, we observed in our study a non¬ 
monotonic dependence of the depletion index on E, with maximum depletion occurring 
for T « 2. In the experiments of Rusconi et al. (2014), a similar non-monotonic trend was 
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reported, with the strongest depletion occurring in the range of 7 « 2.5 - 10 s“^. From 
their data, we estimate Pcf « 5 - 20 and PCg « 0.125, from which we find P « 0.8 - 1.6 
in reasonable agreement with our numerical results. A simple analytical model based on 
a Fokker-Planck equation was also introduced in their paper, though only limited results 
were obtained in the low-Pe/ limit. 

Since the experiments of Rusconi et al. (2014), the existence of centerline depletion in 
strong flows was also confirmed in the numerical simulations of Chilukuri et al. (2014), 
which provided additional insight into the shape of the depletion layer and its scaling 
with flow strength. By fitting the dip in concentration at the centerline with a parabola, 
they were able to extract the profile curvature from their simulation data, and showed 
that it collapses onto a master curve when plotted vs 7 ^, 77 / 214 , in agreement with our 
prediction that the shape of the depletion is controlled by y = Pcg/Pcf = Vs/2jyjH. 
Their also reported similar particle orientations as predicted in figures 6 ( 0 ) and 8 ( 0 ): 
namely, swimmers are preferentially aligned with the flow in the bulk of the channel, 
even though they tend to swim upstream near the walls. Finally, we recall that our 
theoretical scaling for the width of the depletion layer is also in agreement with the 
analytical model of Zdttl & Stark (2012), which is discussed in more detail in Appendix 
D and determines the distance away from the centerline where a deterministic swimmer 
leaving z = 0 with a given orientation fully aligns with the flow, i.e., becomes trapped 
by shear alignment. 


5.3. Concluding remarks 

The favorable agreement of our predictions with both experiments and simulations val¬ 
idates our model and in particular our choice of boundary condition. We reiterate that 
particle-particle and particle-wall hydrodynamic interactions were entirely neglected in 
this work, suggesting that the salient features of confined active suspensions such as wall 
accumulation, upstream swimming, centerline depletion and shear-trapping can all be 
explained in the absence of such interactions. Yet even in dilute suspensions, particle- 
wall hydrodynamic interactions are known play a role (Spagnolie & Lauga 2012) and are 
expected to slightly modify the results described here. Pusher and puller suspensions are 
no longer equivalent when hydrodynamic interactions are included and therefore may 
adopt slightly different distributions, whereas this distinction is irrelevant in the present 
model. As particle density increases, we also expect particle-particle hydrodynamic inter¬ 
actions to become significant, and to destabilize the equilibrium distributions obtained 
in §3 if the concentration is sufficiently high. A preliminary one-dimensional stability 
analysis accounting for flow modihcation by the particles suggests the existence of a 
symmetry-breaking bifurcation above a critical concentration in suspensions of push¬ 
ers, leading to unidirectional flow with net fluid pumping; such an instability was also 
previously predicted using various phenemenological models for active liquid crystals 
(Voituriez et al. 2005; Edwards & Yeomans 2009; Ravnik & Yeomans 2013; Fiirthauer 
et al. 2012; Marenduzzo et al. 2007b). Further increases in concentration may also lead to 
the onset of bacterial turbulence (Marenduzzo et al. 2007a; Gachelin et al. 2014). These 
predictions have yet to be confirmed from a hydrodynamics first-principles perspective 
and may also be investigated computationally using a generalization of the finite-volume 
algorithm presented in Appendix C, or by numerical solution of the approximate equa¬ 
tions for the orientational moments of the distribution function, which were shown to 
be highly accurate in the absence of an external flow. Since the equilibrium states under 
confinement are non-uniform and polarized in the wall-normal direction, the instabilities 
in confined active suspensions could have multifold origins. 

Our study has only focused on the limit of high-aspect-ratio particles whose orienta- 
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tional dynamics are described by equation (2.5). If the aspect ratio of the particles is 
not high, some of the conclusions of this work may change. The distributions in the ab¬ 
sence of flow, including the formation and structure of the wall accumulation layers, are 
not expected to change even in the limit of spherical particles, as confirmed by previous 
simulations of Brownian active spheres (Elgeti & Gompper 2013). However, small-aspect- 
ratio particles will be subject to a weaker alignment with the local shear in an imposed 
flow, which is expected to widen and eventually suppress the centerline depletion layer 
in strong flows. This concept may provide interesting avenues for the sorting of active 
particles by shape in microfluidic devices. 

As a final comment, we recall that a crucial ingredient of our analysis is the presence 
of translational diffusion in the dynamics of the swimmers, which acts to balance the 
swimming flux at the boundaries and leads to diffuse accumulation layers. In the limit 
of strong propulsion or weak diffusion (A —> 0), we saw that accumulation is enhanced, 
and we expect the formation of concentration singularities at the walls in the strict limit 
of dt = 0. This limit is not easily addressed in the context of our theory, though a very 
recent attempt at describing accumulation in this case was proposed by Elgeti & Gompper 
(2015). The development of a more detailed framework in the absence of diffusion may 
prove particularly relevant for describing the accumulation of fast-swimming bacteria 
undergoing run-and-tumble dynamics, notably in applications involving the interaction 
of bacterial suspensions with suspended passive objects (Sokolov et al. 2010; Di Leonardo 
et al. 2010; Koumakis et al. 2013; Kaiser et al. 2014). 

The authors thank John Brady, Anke Lindner, Eric Clement, Roman Stocker, Roberto 
Rusconi and Jeffrey Guasto for useful conversations on this problem. D.S. gratefully 
acknowledges funding from NSF CAREER Grant No. CBET-1151590. 


Appendix A. Comparison between the no-flux and reflection 
boundary conditions 

In this Appendix, we compare the no-flux boundary condition of equation (2.8), which 
is central to our model, to the reflection boundary condition used in previous works 
(Ezhilan et al. 2012; Bearon et al. 2011). The reflection boundary condition ensures that 

= (Al) 


at the channel walls, where 9 and </> are defined in Figure 1. Calculating the first three 
orientational moments of equation (A 1) yields the following conditions to be enforced at 
z = ±1: 




= 0 


dm,. 


dz 


= 0 , 


dz 
dD 


= 0 , 


yy 


dz 


= 0 , 


Dy^ = 0 . 


(A 2) 
(A3) 
(A 4) 


While equations (A2)-(A4) are easily shown to imply that the no-flux conditions (2.25)- 
(2.27) on c, niy^Dyy, D^z are also satisfied, they are much more stringent conditions, with 
a significant impact on the distribution of particles near the wall. 

First, in the absence of flow, we see that equations (3.4)-(3.6) now need to be solved 
subject to boundary conditions (A 2)-(A 4) at z = ±1. The uniform and isotropic solution 
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with = 1 and = 0 satisfies this system exactly. In other words, the 

condition of A 1, by enforcing a zero concentration gradient and wall-normal polarization 
at the walls, is unable to capture the concentration/polarization boundary layer which 
is one of the key results predicted by the no-flux boundary condition and is a ubiquitous 
feature of experiments and particle models. 

The impact of condition (Al) on distributions under flow can be understood in the 
low Pcf limit by modifying the derivation of §4.1. Since = 0, the right-hand term 
in equation (4.7) now vanishes. Equation (4.7)-(4.8) are then rewritten as 


Pe 
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at z = ±1. 


(A 5) 
(A 6) 


(A 7) 


Taking a cross-sectional average of equation (A 5) subject to equation (A 7) shows that 
friy^^ = 0. Therefore, the mean upstream velocity in the channel is exactly zero if the 
reflection boundary condition is enforced. The condition also imposes a zero streamwise 
nematic alignment (DyV = 0) at the walls, which is not physical when a fluid flow 
satisfying the no-slip boundary condition is imposed. A closer look at equations (A6)- 
(A 7) also reveals that the system is in fact ill-posed in the limit of Pcg —>■ 0. For finite 
values of Pcg, a numerical solution shows that the reflection boundary condition severely 
underpredicts the near-wall upstream polarization shown in figure (6). Finally, we note 
that the analysis presented in §4.2 in the strong-flow limit (and hence the scalings for 
the depletion boundary layer thickness and rationalization of the non-monoticity of the 
depletion index with Pcf) describe the dynamics in the bulk of the channel and is not 
affected by the boundary condition imposed. 


Appendix B. Effect of steric exclusion 

The analysis of this paper entirely neglected the finite size of the active particles and in 
particular did not account for steric exclusion with the boundaries, which is expected to 
modify the distributions near the walls as observed experimentally (Takagi et al. 2014). 
As previously shown in the case of passive rods (Nitsche & Brenner 1990; Schiek & 
Shaqfeh 1995; Krochak et al. 2010), excluded volume interactions can be incorporated 
by means of a more complex boundary condition. One must first realize that steric 
exclusion prohibits those configurations near either of the two walls that lead to overlap of 
a section of a particle with the wall. The boundaries between such allowed and prohibited 
configurations define two hypersurfaces in the three-dimensional (z, 6, (f) space of particle 
configurations: 


z = 1 — L* |cos0| (top hypersurface), 
z = —1 -I- A* |cosd| (bottom hypersurface), 


(Bl) 

(B2) 


where L* = LI2PI is the ratio of the particle length to the channel width. At any 
position z inside the channel, this restricts the allowable range of 9 to an interval of the 
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z 


Figure 14. (Color online) Effect of steric exclusion on the steady concentration profile in the 
absence of flow and for Pcs = 0.25. The plot compares numerical results for three different 
valnes of L* = L/2H to the case where steric exclusion is neglected (L* —>■ 0). 
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and consequently, any integral with respect to p of a field variable A(z,p) must be 
restricted to these configurations: 


p ^ 27 t ^ 62 ( 2 ) 

/ A{z,p)dp= / / A{z,p)sm.Q dO dcf) 
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To ensure that prohibited conhgurations are never realized, the boundary condition 
(2.7) must be replaced by a more general no-flux condition on the hypersurfaces defined 
in equations (B 1)-(B 2). Introduce the generalized flux vector J as 


with 
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Denoting by n(z,9) the normal unit vector on one of the two hypersurfaces, the gener¬ 
alized no-flux condition is simply expressed as 


n(z,9) •J{z,p,'I') = 0, 


(BIO) 
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which, upon calculation of the normal n, leads to the two conditions: 

L*sinejg = 0 at z=l-L*\cos9\, (B11) 

J^±L*sin6lJe = 0 at z =-1 + L* |cos6»|. (B12) 

In each case, the upper sign is used when 9 G [0,7t/ 2] and the lower one when 9 G 
[7t/2,7t]. Numerical solution of the conservation equation (2.6) subject to the boundary 
conditions (B 11)-(B 12) can be done using finite volumes as described in Appendix B. 
Typical results for the concentration profile c(z) in the absence of flow are shown in 
figure 14 for different values of L* and compared to the solution obtained previously 
using the boundary condition (2.7), which corresponds to the limit of L* —>■ 0. When 
steric exclusion is accounted for, a depletion layer is observed close to the walls whose 
thickness is of the order of L*. Steric exclusion leads to a decrease in concentration in 
the near wall region because it suppresses the orientations aligned towards the wall and 
hence the wall normal polarization. Under stronger confinement (higher L*), this leads 
to a concentration peak at the edge of the depletion layer due to wall accumulation, and 
this peak increases in magnitude and shifts closer to the wall as L* decreases. For very 
small values of L*, the concentration profile approaches the profile obtained by neglecting 
steric effects, and steric exclusion can be safely neglected outside of the depletion layer 
itself whenever L* < 0.01. This is indeed the appropriate regime in most microfluidic 
experiments with bacterial suspensions, which justifies the use of the simpler boundary 
condition (2.7) in the work presented here. 


Appendix C. Finite-volume numerical algorithm 

In this Appendix, we describe the algorithm used for the numerical solution of equation 
(2.13) for the distribution function. The method is based on a finite-volume discretiza¬ 
tion of the Smoluchowski equation (Ferziger & Peric 2002), which has the advantage of 
satisfying conservation locally to machine precision while also allowing for an easy im¬ 
plementation of no-flux boundary conditions such as (2.7) or (B 11)-(B 12). To avoid the 
cost of large matrix inversions, we solve the time-dependent Smoluchowski equation to 
steady state using an explicit scheme. In conservative form, the governing equation can 
be written as 

dT 

_ + V,7.J = 0, (Cl) 

where J is the generalized flux vector defined in equations (B6)-(B9), and Vj is the 
gradient operator in the three-dimensional (z, 0, 4>) space of particle configurations: 
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We note that T{z, 9, (f) is defined on a hypervolume obtained by extruding the unit sphere 
in the z dimension. This computational domain is discretized into finite volumes using a 
uniform grid with respect to (z, r, (f), where r = cos 9. The nodal points (z*, r^, (ff) where 
is evaluated are located at the centers of each volume and have coordinates 




2i-l 

N. 

- 1 

for 1 = 1,. 


(C3) 

2j-l 

Nr 

- 1 

for j = l,. 

■ •) AC, 

(C4) 

2n{k — 
Ns 

i) 

for k = 1, 

...,N^, 

(C5) 
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Figure 15. Typical finite volume in three-dimensional {z,0,(j>) space, centered around an ar¬ 
bitrary nodal point with indices The uniform discretization with respect to {z,r,4)) 

ensures that all such computational cells have equal volume AV — AzArAcf). 


where Nz, Nr, and are the total numbers of points in each direction. We also define 
the grid spacing in each direction as 




Alt) 


2n 


(C6) 


The advantage of this discretization (compared to a uniform grid with respect to 0) 
is that it divides the sphere of orientations into elements of equal area, which reduces 
restrictions on the time step arising from the rotational flux. 

A typical finite volume centered around node {i,j,k) is illustrated in figure 15. It is 
delimited by eight grid points denoted A through H, with indices (i±,j±, k±) where we 
have introduced the notations i± = f ±0.5, j± = j±0.5 and k± = k±0.5. The cell edges 
have lengths 


AB = DC = EF = EdG = Ale = cos — cos ^(r'^+). 


AD = EH ^ All = 


27Tsind'l- 


BC = EG = Alt = 


2n sin 0-1+ 


AE = BE = DH = CG = Az. 


(C7) 

(C8) 

(C9) 


In figure 15, faces ABCD and EFGH have unit normal z and surface area ArAij). 
Similarly, faces ADHE and BCGF have unit normal 9 and areas AzAl^jt and AzAl'^, 
respectively, whereas faces ABFE and DCGH have unit normal cj) and area AzAlg. 
The volume of the computational cell is AV = AzArAij). 

In order to satisfy conservation of the distribution function exactly in each finite vol- 








Transport of a dilute active suspension 
ume, we first integrate equation (Cl) over computational cell V(i,j,k): 

dT 


—— h V j • J dz dr d(/) = 0. 
ot ) 

V(,i,j,k) 

After applying the divergence theorem to the second term, this can be recast as 
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d 

0 = — I 11 T dz dr d(l> + 

ot ' ' ' 


Jz dr d(j) — 


Jz dr d(j) 


Vd,J,k) 


ABCD 


EFGH 


+ / / Jedzdcj) — / / Jgdzdcj) 


(Cll) 


ADHE 


BCGF 


j j J(p dz dr — JJ J(j, dz dr. 


+ 

ABFE 


DCGH 


Volume and surface integrals in equation (C 11) are approximated to second-order using 
a midpoint rule. After division by AV, this leads to the discretized equation: 


0 = 1 “ \ k) — Jz{i—,j, k)] 


1 
Ar 
Ale 


+ TC; [Je{'i,j+,k)s\TLe^+ - Je{i,j-,k)sine^-] 
[J4,{i,3, k+) - J4>{i,j, k-)]. 


(C12) 


ArA(j) 

In order to integrate this equation, we must first obtain approximate expressions for the 
fluxes at the centers of the six volume faces. This is done using linear interpolation for 
terms involving T, and centered finite differences for terms involving derivatives of T. In 
the z and <j) directions, this gives 


Jzii+,j,k) « PCs cos 9^ 


q,i+l,j,k _|_ yjyi,j,k 


- 2APei 




Az 


J^{i,j,k+) « i 


PefS{z^) cos 9^ cos (j>^ 


q,i,j,k+l _|_ q/i,j,k 


sin 93 


IpiJA+l _ ^p^,j,k\^-^ 

Af) 


(C13) 

(C14) 


with similar expressions for Jz{i-,j,k) and J^{i,j,k-). The approximation of Jg is 
slightly more involved due to the non-uniformity of the mesh with respect to 9. Deriva¬ 
tives with respect to 9 are calculated using symmetric central finite differences in terms 
of r after application of the chain rule, and linear interpolation is used with respect to 
the 9 variable, leading to the approximation 


Jg{i,j+,k) « i {Pe/S'(z*)cos0J+ coscj)'^ [A^+ifbi+i.fc 


sin 9^+ 


gpi,j + l,k _ gpi,j-l,k 

Ar 


(C15) 


with a similar expression for Jg{i,j-, k). The interpolation weight A'’+ is given by 

A^+ = 


cos ^(r3 + — cos ^{r3) 


cos~^ {rj + Ar) — cos~^ {r3)' 


(C16) 
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When integrating equation (C 12) in time, care must be taken when dealing with cells 
adjacent to the poles of the unit sphere (j = 1 and N^), as these cells are missing one 
face. For instance, cells with j = 1 are such that A = D and E = H in the diagram 
of figure 15, so that face AD HE is missing and the corresponding flux should not be 
included in the discretized equation. 

Boundary conditions also need to be specified to proceed with the time integration. 
Periodic boundary conditions are used in the (j) direction, yielding: 

J^il/2,j,k) = 1/2, j,k) and + 1/2, j, k) = J,p{3/2, j, k). (C 17) 

Treatment of the boundaries in the 9 and 2 ; directions differs depending on whether steric 
exclusion with the walls is included or not. 


C.l. Without steric exclusion 

When steric exclusion is not included and the simple boundary condition of equation 
(2.7) is used, 6 varies over its full range [0, tt]. However, no boundary condition is needed 
along 6 as the boundary cells with j = 1 and Nr are missing one face as explained above, 
which eliminates the need to specify Jg{i,1/2,k) and J,f,{i,Nr + 1/2,k). Along the z 
direction, the boundary condition is simply the no-flux condition (2.7), which translates 
into 

JS,j, 1/2) = Mi,j,N, + 1/2) = 0. (C 18) 


C.2. With steric exclusion 

The situation is more complex when steric exclusion is accounted for, as the boundary 
conditions needs to be enforced on the hypersurfaces defined in equations (B 1)-(B 2). It 
is convenient in this case to choose and Nr such that 

Z\z = L*Ar or N^ = —. (C 19) 

L* 


Indeed this ensures that the hypersurfaces fall onto grid points and eliminates the need 
for further interpolation. However, if L* is small, this implies that a significantly finer res¬ 
olution is needed along z than along 6. As we discussed in Appendix A, the hypersurfaces 
limit the range of allowable values of 9 to an interval of the form [9i(z),92(z)] C [0,7t] 
for particles located near the walls. After discretization of the domain and choosing N^ 
and Nr to satisfy condition (C 19), we find that for any nodal point with coordinate z*, 
there is a finite range gP 2 (*)] of allowable values of 9^, with 


Ji(*) 


h{i) 


Nr 


^ + 1 -* 

if z < — 1 -1- 

2 


Nr 

if z > 1 — L*, 

1 

otherwise, 

Nr 

if z < -I 

2 


Nr 


-i+N. + l 

— i ifz>I — 

Nr 

otherwise. 


(C20) 


(C21) 


Interior nodal points such that j G [ji(*) + l,j 2 {i) — 1] are such that full cuboidal finite 
volumes in (z, r, (f>) can be constructed around them, and therefore do not require any 
special boundary treatment. Boundary nodal points such that j = ji{i) or j 2 ii), however. 
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are contained inside prisms whose hypotenuses coincide with the hypersurfaces. These 
hnite volumes can be treated in the same way as interior control volumes by prescribing 
zero-flux contributions from surfaces lying outside of the domain, by multiplying the 
volume AV by 0.5, and by adjusting the surface area of faces ABFE and DCGH to a 
reduced triangular area given by 


AA = 


AzAlg 

2 


G Alg — 2 sin 



cos 


03 + _|. 03 - 


(C22) 


Appendix D. Active particle trajectories and shear trapping 

In this Appendix, we rationalize the linear dependence of the depletion layer thickness 
5d upon Pcs/Pef by deriving the trajectory of a deterministic swimmer whose dynamics 
result from self-propulsion and shear rotation via Jeffery’s equation. A similar derivation 
was previously presented by Zottl & Stark (2012, 2013). In dimensional variables, the 
equations of motion of the swimmer are written 

z{t) = Vs cos 9 (t), (D 1) 

p{t) = {l-pp).{CE + W).p. (D2) 


Here, C is a shape parameter, with C « 1 for a slender particle as we have assumed in 
the rest of the paper. The two second-order tensors E and W are the rate-of-strain and 
vorticity tensors of the imposed flow, respectively: 

E = + zy), W = ^z{t) {yz - zy) . (D 3) 

Parameterizing the orientation vector as p = (sinycost/), sinysin( 5 i, cos0), we can use 
equation (D 2) to obtain expressions for the time rate of change of the polar and azimuthal 
angles of the swimmer as 


m 

^{t) 


'^z{t)sm(j){t) [(^ -I- 1) cos^ 6{t) — (C — 1) sin^ d{t)] , 

cos9{t)cos(l){t) 

+ sm9{t) ■ 


(D4) 

(D5) 


Equations (Dl), (D4) and (D5) form a closed system of coupled ordinary differential 
equations that can be solved for the swimmer dynamics. 

Any swimmer that is not perfectly aligned with the walls (cos 6* ^ 0) will tend to 
migrate towards one of the boundaries due to self-propulsion, while shear rotation tends 
to align it along the flow direction causing it to get trapped. Recalling the definition of 
X as the ratio of the time scale for shear rotation to the time it takes for a swimmer to 
cross the channel. 


Pef ’ 


(D6) 


we expect two different regimes. When y ^ 1, any swimmer released from the centerline 
with initial orientation (9q, (fo) will reach one of the walls before becoming trapped. On 
the other hand, when y <C 1, we expect there to exist a position Ztrap{9oj ^o) inside the 
channel where the swimmer gets trapped due to shear alignment. This indeed corresponds 
to the regime discussed in §4.2, where depletion from the centerline and shear-trapping 
were predicted to occur for Pcg <C 1 and Pe/ ^ 1. 

To derive a quantitative estimate for Ztrap, we calculate the value of 2 : at which 9 first 
reaches ±7 t/ 2. We first consider the case of a particle with initial position zq = 0 and 
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orientation defined by 9q G [0,7t/ 2), = 37t/2. For this specific initial configuration, 

(/)(0) = 0 which implies ^(t) = 37 t/ 2 for all times. The motion is two-dimensional in this 
case, and the dynamics is governed by the two coupled ordinary differential equations 


z{t) = Vs cos 9 (t), (D7) 

^(0 = [(C + 1) cos^ - (C - 1) sin^ 0{t)] . (D 8) 


An equation for the swimmer trajectory can then be obtained by taking the ratio of (D 7) 
and (D 8 ): 


d9 z (C-b 1) — 2 — 1) sin^ 6 l(t) 

dz H 2xcos0 

This can be integrated from {z,9) = ( 0 , 6 * 0 ) to {ztrap,T^I‘2'), yielding 


(D9) 


^ Ztrap{9o) '^^ _ 


2C(C+1) 


tanh 


^ — tanh ^ 


C + 1 


I 2C 
C + 1 


sm 9o 


(DIO) 


For a typical swimmer of aspect ratio 10, we estimate C « 0.98. Taking the initial con¬ 
figuration to be 6*0 = 0, equation (D 10) simplifies to Ztrap/H ss ^3% « l.li^Pes/Pcf. 
This estimate is consistent with the high-Pe/ scaling analysis of §4.2, as well as with the 
numerical results of §4.3 where we found 6d ~ 2.404 Pcs/Pcf. 

The more general case of an arbitrary initial orientation (9Q,(j)Q) can also be solved 
analytically. Combining equations (D4) and (D 5) to eliminate z{t), we find after inte¬ 
gration: 


cos (j) = cos ( 


(C + 1) cosec^ 6 * — 2C, 


(C + 1) cosec26*o — 2Z 
Now, using equations (D 1) and (D4), we get 


(Dll) 


/ Ztrap{4'0^^0^ 
\ H 



cos 6 * 

(C + l-2Csin2 9) yr — COS^ (j) 


d 6 », 


(D12) 


where sin^’ is known in terms of 9 using (D 11). This expression confirms the scaling of 
Ztrap with yZC) and it can in fact be shown that Ztrap in equation (D 12) has an upper 
bound given by the previous estimate (D 10). 
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